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SUMMARY 

A new crack element has been developed and incorporated into COSMIC/ 
NASTRAN. The element is considered linear, isotropic, and homogeneous. Mode 
I and II stress intensity factors are automatically calculated. Comparisons 
to theoretical p^ane strain solutions for several geometries are presented and 
demonstrate the accuracy of the developed element. Extensions of the element 
'n three dimensions, anisotropic material, and plastic analysis are discussed. 


INTRODUCTION 

Crack or singular elements have been developed for finite element analysis 
for almost as long as finite element codes have been available* These elements 
usually are classified as either hybrid or singular element formulations. Many 
of the elements developed suffered from either lack of accuracy, generality, or 
consistency. Barsoum (Ref. l) points out shortcomings of several different 
elements. These shortcomings include inability to mod^l rigid body or constant 
strain modes, inability to include thermal or body force effects, and lack of 
compatibility with other elements. 

The elements developed by Barsoum (Ref. 1) and Henshell (Ref. 2) rectified 
many of the p oblems described above; however, these elements were limited to 
displacement of the form r^'^. Consequently, they could only model strain 
singularities of the form r _1 ' Recently, Stern (Ref. 3) and more recently, 
Hughes and Aikin (Ref. 4) have introduced families of consistent, conforming 
elements which allow displacements of the form r Y . While the Stem element 
appears to have the restriction that 0 < y < 1 , the element of Hughes and 
Aikin is valid for all y > 0. The element described herein is based upon 
shape functions suggested by Hughes and Aikin (Ref. 4). 

The element presented here possesses the required rigid body and constant 
strain modes. It properly models thermal, body force, and pressure loading 
conditions- Additionally, it is compatible with standard linear or quadratic 
isoparametric elements and can possess either 5 or 6 nodes. Finally, it can 
<s used as a nonsingular element with a variable number of nodes. 
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ELEMENT FORMULATION 


The following derivation follows Hughes and Aikin (Ref. 4). Referring to 
Figure 1, the standard bilinear shape functions are used for nodes 1 through 
4: 


N^r.s) * (i-r)(l-s) 

Not 1 "’ 3 ) = I ( 1_ s) 

(D 

N-(r,3) = rs 
R 4 (r,s) = (1 -r)s 


The shape functions for nodes 5-8 are chosen as: 


N 5 (r,s) = (l-s)P(r.T) 

Ng(r,s) = rP(s,2) 

( 2 ) 

N 7 (r,s) * sP(r, y) 

N 6 (r,s) = (l-r)P(s,2) 


where 


F(x,y) = 2( x - 


x* - jj 1J2]\ } 

1 - 2(l/2) Y 


(5) 


It can be easily shown that the shape functions for nodes 5-8 reduce to 
the standard quadratic serendipity element when y of Equation (5) is set equal 
to 2. It can also be seen that the shape function for nodes 5-8 satisfies the 
interpolation property at all nodes of the element. That is: 


N.(r.) = 6. . and N (s .) = 6 . 
i 3 ij i j ij 


where r. and are values of r and s at node j and 6. . is the Kronecker 

j j i i 

delta. However, the shape functions associated with nodes 1-4 do not satisfy 
the interpolation property at nodes 5-8. Following the standard technique 
(Ref. 4), tne shape functions fo; nodes 1-4 are modified as follows: 
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N 1 ♦ N 1 (r,s) - [N 8 (r,s) 
« 2 ♦ N 2 (r,s) - [N^ (r,s) 
N ? ♦N 3 (r,s) - [N 6 (r,s) 
N 4 ♦N 4 (r,s) - [» 7 (r,s) 


♦ N 5 (r,s)]/2 
+ N 6 (t,s)]/2 
+ N„(r,s) ]/2 
- Kg(r,s)J/2 


where the + reads: "is replaced by." 


( 4 ) 


It can now be oeen that the shape functions for all eight nodes satisfy 
the required interpolation property. Additionally, the sha^e functions are 
capable o** exactly representing the monomials 1 , r,s,r\rs,s ,r T s, and s^r. 

The ^ eaence of 1 ,v and s ensure representation of rigid body and constant 
strain modes. The presence of r Y allows exact representation of displacements 
of the # f^rm r Y . Note that this will result in a line singularity of the 
fora r‘~ upon differentiation. 


In order to represent point singularities, the quadrilateral must be 
degenerated into a triangle. This is done by coalescing nodes 4, 8, and 1 as 
can be .one for standard isoparametric elements (Eef. 5) and as is shown 
schematically in figure 2. Thus, finally, for a point singularity, the shape 
function associated with node i is replaced with: 


N lV 'r,s) ♦ (r.s) + N^r.s) ♦ S 8 (r,s) (5) 

This is easily programmed into the element routine. The generality of 
the above derivation allows use of the same shape functions as the basis of 
three-dimensional elements which possess line singularities. 

In summary, for the 6 node triangle, the shape function associated with 
node 1 is given by Equation (5), the shape functions associated with nodes 2 
and 3 are given by N 2 and N-j cf Equation (4), and the shape functions 
associated with nodes 5 through 7 are given by N^ through of Equation (2). 

Given the shape funcuions for the element, calculation of the stiffness 
aatrtr, thermal load vector, and gravity load vector follows the standard 
procedure as described in Reference 6. These quantities are therefore given 
as: 
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K® = J B T D B dV 
v e 

= / B T E a AT dV 

V e 

e r T 

F^ - N‘ bt iV 
~b * *'* *■* 


( 6 ) 


where 


B = L S , N = [N I, N 2 I, ...] 


3_ 

3x 


L 


0 


0 

3y 


3 3 

Uj "57 


and I is a 2 x 2 identity matrix. 


See Reference 6 for more details. The actual integration is performed via 
Gaussian quadrature. That is, the integrals are approximated as: 


/ f(x) dx = l V. f (a .) (7) 

j-1 J J 


Due to the formulation, it can be shown that along the s direction, the 
integration order needs to be, at most, 4 in order to exactly integrate the 
element. For an undistorted element, a naxinmm integration order of 5 is nec- 
essary, although often an integration order of 2 yields results just as good. 
Along the r direction, the method for an exact integration has not been ascer- 
tained. Currently, an integration oruer of 4 or 5 seems to suffice along the 
r direction. An exact integration formula analogous to the formula presented 
in Reference 3 or Reference 7 will hopefully be derived in the near future. 

Calculation of the stress intensity factors are performed using the 
equations: 
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( 8 ) 



where the nomenclature is shown in Figure 3* Alternatively, similar equations 
in terms of stresses can be used. It has been found here and noted elsewhere 
(Ref. S) that Equation (8) yields more accurate results than the similar 
equations in terms of stresses. The values of the stress intensity factors 
are calculated at each of the integration points along the r direction and 
extrapolated to r=0 using Lagrangian interpolation. 

Extension of the above formulation to three dimensions is straightforward. 
As discussed in Reference 4, the three-dimensional shape functions are simply 
products of the two-dimensional shape functions in r and s with the desired 
one-dimensional shape function in t. For example, the shape functions for the 
three-dimensional element of Figure 4 are given by: 


N. (r,s,t) 


tNj (r,s) 


N 6 (r,s,t) 
N 7 (r,s,t) 

« 


tNgtr.s) 

(l -t)Nj (r,s) 
0-t)N 6 (r,s) 


The above element formulation may be extended to anisotropic materials by 
using the appropriate anisotropic material matrix D in Equation (6). When 
appropriate* D could be different at each integration point. 

In order to incorporate plasticity effects, it is suggested that standard 
techniques currently used for plasticity in regular elements could also be 
applied to the present element. That is, after each load increment, each 
integration point in the element would be checked to see if it has gone plastic 
or not. If plasticity has occurred, then an algorithm such as radial-return 
(Ref. 9) would be used to bring the stress back to the yield surface. The 
element's internal forces, used to calculate the out-of-balance load vector 
would be given in standard form as: 
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e 



0 d V 


Additionally, the order of the strain singularity would have to be updated, 
depending upon the hardening properties of the material (Ref. 10). 


IMPLEMENTATION IN NASTRAN 

Implementation of the creek element into NASTRAN was performed via the 
dummy element CDUM1 . This procedure is covered in Reference 11. The present 
element was modeled after the QDMEM1 element routines, due to their similarity. 

The first step was to create a subroutine KDUM1 which generates th'? 
stiffness and mass matrices. The mass matrix may be either consistent or 
lumped. When the mass matrix is used for calculation of gravity loads, the 
consistent mass matrix should be specified. This subroutine is eventually 
linked to NASTRAN LINK 8. 

For computation of thermal loads, the subroutine EDTL must be modified to 
make a call to SSGETD before calling the routine DUM1. The dunuqy coding in 
routine L’JMI is then modified to calculate the thermal load vector based on 
the connecting grid point temperatures. Optionally, the element centroidal 
temperature could h^-e been used, although this is generally not recommended 
since temperature gradients near t v e crack could not be accurately represented 
in this way. After modifying EDTL and DUM1 , they must be linked to NASTRAN 
LINK 5. 

Finally, the dummy coding for the SLUM11 and SDUK12 routines must be 
modified so that they perform the required operations. SDUM11 performs the 
preliminary geometry calculations and creates the S matrix which relates 
element stresses (including stress intensity factors) to the element’s grid 
point displacements. SDUM12 then uses the S matrix, grid point displacements, 
and temperatures to compute the centroidal stresses and stress intensity 
factors and writes them to the output file. After modifying the SDUK11 and 
SLUM -i 2 coding, it is linked to NASTRAN LINK 13. This completes the 
implementation into NASTRAN. 


NUMERICAL RESULTS 

In order to assess the accuracy of the present element, four different 
crack geometries/loading conditions with known solutions were analyzed. 

Figure 5 shows the different geometries analyzed. Figure 6 presents four 
different mesh sizes which were used to analyze the first three crack geome- 
tries. Figure 7 shows the boundary conditions used. For the edge crack with 
a point load, Figure 7 is modified so that the load is applied at the edge of 
the crack. Table 1 presents the errors associated with both the crack opening 
displacement (COD) and the mode I stress intensity factor Kj. As can be seen, 
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the COD is less sensitive to the mesh size, while the Kj values appear to be 
converging to their exact solutions. However, the edge crack with a point 
load solution appears to overshoot the exact by about 5$. It should be men- 
tioned that the "exact” solution for the edge crack specimen with a point load 
is considered to be accurate to within 2 % . The other exact solutions were 
considered to have accuracies better than These exact solutions were 

obtained from Reference 12. 

Figure 8 presents a model of a central crack in a finite plate. To ascer- 
tain the accuracy of the element^ mode II stress intensity factor, Kjj, the 
model of Figure 8(a) was used. The results for both Kr and Kjj are presented 
in Table 2. As can be seen, the K n is within about 4% cf the exact solution. 

In summary, both COD and stress intensity factors appear to be accurately 
represented even for relatively course meshes. The accuracies obtained are 
well within the accuracies required by typical engineering calculations. This 
is due to the fact that the scatter alone, in the Kj values during a typical 
test, may be 1 0^ . 


CONCLUSIONS AND FUTURE RESEARCH 

An element formulation has been presented that accurately models singu- 
larities. The form of the singularity is general and the two-dimensional 
element developed may be easily extended to three dimensions. Additionally, 
the element may be used as a standard, variable number of nodes quadratic 
element. The element has been incorporated into NASTRAN and compared to 
several known solutions. 

Future research will include more tests of the element against known 
exact solutions. Additionally, an exact integration rule is desirable, and 
work to develop this will be performed. The element formulation will then be 
extended to three dimensions and the code will be incorporated into NASTRAN. 
Finally, extensions to include anisotropic materials and plasticity are 
possible and should be studied further. 
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TABLE 1 ERRORS IN COD AND Kj AS A FUNCTION OF MESH 




37 Grid Mesh 

77 Grid Mesh 

86 Grid Mesh 

95 Grid Mesh 

CENTRAL CRACK 

COD Error (%) 

- 7.73 

- 2.46 

- 2.82 

- 2.92 


Kj Error (%) 

11.63 

2.91 

1.69 

1.22 

EDGE CRACK WITH 
UNIFORM STRESS 

COD Error (%) 
Kj Error (%) 

- 5.24 

- 6.26 

- 2.51 

- 5.18 

- 2.60 
- 3.19 

- 2.64 

- 2.05 

EDGE CRACK WITH 
POINT LOAD 

COD Error (%) 
Kj Error (%) 

-26.11 

- 9.38 

- 0.40 

5.26 


TABLE 2 ERRORS IN Kj and K n FOR 234 GRID MESH 




Error (%) 

CENTRAL CRACK IN 
SHEAR, FINITE PLATE 

K I 

- 1.37 


K II 

4.23 


























s r=l 



1 5 2 


Figure 1 Nomenclature for Eight Node Isoparametric Element 
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Figure 2 Degeneration of the Eight Node Element to a Six Node 
Triangular Element 
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crack 

Figure 3 Nomenclature for Crack Geometry 



Figure 4 A Possible Three-Bimensional Generalization of the 
Two-Dimensional Element of Figure 2 
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Figure 6 Different Mesh Sizes Analyzed 

(a) 37 grid mesh 

(b) 77 grid mesh 

(c) 86 grid mesh 

(d) 95 grid mesh 



(a) (b) 


Figure 7 Boundary Conditions for Edge Crack and Central Crack Specimens 

(a) edge crack 

(b) central crack 


1 h 

(a) loading condition for Kjj 


calculation 



(b) loading condition for Kj calculation 
Figure 8 Model of Central Crack in Finite Plate 
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